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manifold using implicit DEC scheme 
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Abstract 

Maxwell's equations can be solved numerically in space manifold 
and the time by discrete exterior calculus as a kind of lattice gauge 
theory. Since the stable conditions of this method is very severe re- 
striction, we combine the implicit scheme of time variable and discrete 
exterior calculus to derive an unconditional stable scheme. It is an 
generation of implicit Yee-like scheme, since it can be implemented in 
space manifold directly. The analysis of its unconditional stability and 
error is also accomplished. 
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1 Introduction 

The Yee scheme is a commonly employed efficient approach to solve Maxwell's 
equations numerically and so to model wave propagation problems in the 
time domain [1]. Although it is not a high order method, it is still preferred 
for many applications because it preserves important structural features of 
Maxwell's equations [2-5]. Bossavit et al present the Yee-like scheme to ex- 
tend the Yee scheme to unstructured grids. This scheme combines the best 
attributes of the finite element method (unstructured grids) and the Yee 
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scheme (preserving geometric structure) [6,7]. Stern et al [8] generalize the 
Yee-like scheme to unstructured grids not just in space, but in 4-dimensional 
spacetime by discrete exterior calculus(DEC) [9-20]. This relaxes the need 
to take uniform time steps. Based on these results, we generalize the Yee- 
like scheme to space manifold and the time, which is a kind of lattice gauge 
theory [21]. 

The stable conditions for Yee-like scheme and its generation are very 
severe restriction, and imply that very many time steps will be necessary 
to follow the solution over a reasonably large time interval. In computer 
simulations of physical processes, explicit and implicit methods are often 
used. For some problems, it takes much less computational time to use 
the implicit method with larger time steps, even taking into account that 
one needs to solve equations at each step. Based on these considerations, 
Keraen et al present the unstable conditional implicit Yee-like scheme [22]. 
In this paper, we show that the implicit scheme and DEC can be united to 
find an unconditional stable scheme (IDEC) for solving Maxwell's equations 
on space manifold and the time. The analysis of IDEC's unconditional 
stability and error is also accomplished. This scheme reduces to the implicit 
Yee scheme, if choosing rectangular mesh for flat space, and reduces to the 
scheme presented by Keraen et al, if choosing tetrahedral mesh for flat space. 

2 Preliminaries 

A discrete differential /c-form, k £ Z, is the evaluation of the differential 
fc-form on all /c-simplices. Dual forms, i.e., forms evaluated on the dual cell. 
Suppose each simplex contains its circumcenter. The circumcentric dual cell 
D{gq) of simplex gq is 

D(g ) := [J Int(c(c7 )c((7i) ■ ■ ■ c(oy)), 

where G{ is all the simplices which contains Go,..., G{-\, and c(<7j) is the cir- 
cumcenter of <7j. In DEC, the exterior derivative d is approximated as the 
transpose of the incidence matrix of Ai-cells on k + 1-cells, and the approx- 
imated Hodge Star * scales the cells by the volumes of the corresponding 
dual and primal cells. 

The 2D or 3D space manifold can be approximated by triangles or tetra- 
hedrons, and the time by line segments. Discrete connection 1— form or 
gauge field A assigns to each element in the set of edges E an element of 
the gauge group E: 

A : E -> E. 
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Discrete curvature 2— form is the discrete exterior derivative of the discrete 
connection 1— form 

F = dA : P M. 

The value of F on each element in the set of triangular P is the coefficient 
of Holonomy group of this face. The 2— form F automatically satisfies the 
discrete Bianchi identity 

dF = 0. (1) 
For source case, we need discrete current 1— form J. Let A = ^ M and the 

E 

Lagrangian functional be 

L(A,J) = -\{dA,dA) + (A, J), 

where 

{dA,dA) := (A) lx | E |(d)| jB | x | F |(*)| F | x | F |(d)| 2 F|x|jB| (^)j J E|xl 
(A, J) := (-4)ix|E|(*)|£'|x|s|(^ T )|E|xi- 

The Hamilton's principle of stationary action states that this variation must 
equal zero for any vary of Ai, implying 

d T * dA = *J (2) 

Since (d T ) 2 = 0, the discrete continuity equation can express as: 

d T * J = 0. (3) 

The Eqs.(l-3) are called discrete Maxwell's equations, which are invariant 
under gauge transformations A — > A + df for any 0- forms. 



3 IDEC for Maxwell's Equations 

If allowing for the possibility of magnetic charges and current discrete 3— form 
J, the symmetric discrete Maxwell's equations can be written as 

dF = J d T *F = *J, (4) 

with discrete continuity equations or integrability conditions 

dJ = d T * J = 0. (5) 
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Implicit scheme for TE wave 

The discrete current forms, discrete curvature 2— form, and their dual can 
be written as 

J n = (pi, -J™ h a dt) *j n = (- * ( Pe dty, *j e n+l ) 

F n = E n+1 A dt + B n *F n = H n+1 Adt - D n , 
where n and n + ^ denote the coordinate of the time, E = ^ £^e* (electric 

E 

field) is the discrete 1— form on space, B = ^BiP 1 (magnetic field) is the 

p 

discrete 2— form on space, H = J2Hi*P l (magnetizing field) is the dual of B 

p 

on space, D = J2 A * e% (electric displacement field) is the dual of E on space, 

E 

p e dt (charge density) is the discrete 1— form on time, J e = Y1 Jei& 1 (electric 

E 

current density) is the discrete 1— form on space. p m = Yl PmiT 1 (magnetic 

Tet 

charges) is the discrete 3-form on space, J m = JmiP 1 (current) is the 

p 

discrete 2— form on space. 

The symmetric discrete Maxwell's equations (4) and (5) can be written 

as 

d B n — n n 

"s-c — P m 



d s E n+1 Adt = -d t B n - Jra~ 2 A dt 

djD n = *(p e dt) n 
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(6) 

djH n+1 Adt = djD n + *jr^ 2 , 
where d s , df are the restriction of d, <i T on space, and 

nn+l _ Dn nn+l _ rya 

d t B n := A dt djD n := A dt. 

At 1 At 

Proposition 3.1 // the initial condition satisfies the first and third equa- 
tions in Eqs.(6), the solution of the second and fourth equations in Eqs.(6) 
automatically satisfy Eqs.(6). 

Proof. Because the dimension of spacetime is 3 + 1 or 2 + 1, therefore 

df * (p e dt) n = df * Je +2 ~ = d sP n m = dt A dt = 0, 
and the continuity equations can be reduced to 

df * ( Pe dt) n -df* JT^ =0 - d s J^ A dt + d t p n m = 0. 
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So we have 



dJd T s D n - df * ( Pe dt) n 
d t d s B n - d t pl 



-df * ( Pe dt) n - d T s (d T s H n+1 Adt- *j f 



n+i 







n 

J rn 



d t p n m + d s {d s E n+1 Adt+J"m*Adt) 



0. 



□ 

Now we show the implicit scheme (6) on the 2D discrete space manifold 
and the time. Take Fig.l as an example for a part of 2D space mesh, in 
which e±,..., e§ are edges, Pi, P2 are triangles, *e\ is the dual of e\. The 
second and fourth equations in Eqs.(6) based on Fig.l are 



n+l 



At 

B n+1 



r J el 



H n+l 



* ed 



At 



E n+1 



|ei| +E2 +1 \e 2 \ +E% '-|e 3 | 



?n+l I 



(7) 



Pi 



where | | denotes the measure of forms and dual. The summation on the 
right is orient, that is to say, inverse the orientation of ej, then multiply — 1 
with Ei. Notice that a significant difference from Yee-like scheme is that 
*6i is a polyline. Eqs.(7) can be implemented on 2D discrete space manifold 
directly (see Fig. 3) so is a generation of implicit Yeedike scheme. 





Figure 1: edge and face with direction 
In the absence of magnetic or dielectric materials, there are relations 

Di = eoEi Bi = /M)Hi, (8) 

where £0 and p,Q are two universal constants, called the permittivity of 
free space and permeability of free space, respectively. With relations (8), 
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Eqs.(7) can be rewritten into an implicit scheme (9) for TE wave. 



1 At +Jel " |*ei| 



rxn+l Tjn , i J7 ,n + 1 U I _i_ T? n +^-\^ I _L P"+lL I 

#1 - ijj n+i i% |ei|+-fe 2 l e 2|+-fe 3 |e 3 | 
MO 7- r J ml 



TE (9) 



Implicit scheme for TM wave 

If writing 

pn = H n+1 Adt _ D n ^pn = _pn+l A df _ gn 

J"" = {-P n e , Je n+I A dt) * J = (- * (p m dtr, *J™ + ^), 

where H = H^e 1 is the discrete 1— form on space, D = ^HiP 1 is the 

E P 

discrete 2— form on space, E = ^ Di * P l is the dual of D on space, B = 

p 

J2Bi*e l is the dual of H on space, p e = J2 Pe(E % is the discrete 3— form on 

E Tet 

space, J e = JeiP 1 is the discrete 2— form on space, p m dt is the discrete 
p 

1— form on time, J m = Y1 JmiE 1 is the discrete 1— form on space, the discrete 

E 

Maxwell's equations can be rewritten as 



i 



d s D n = pl 

d s H n+1 A dt = d t D n + 4 t+5 A dt 

d T s B n = *(p m dt) n { } 

d^E^Adt = -dfB n -*J^. 



Proposition 3.2 If the initial condition satisfies the first and third equa- 
tions in Eqs.(lO), the solution of the second and fourth equations in Eqs.(lO) 
automatically satisfy Eqs.(lO). 

Proof. Because the dimension of spacetime is 3 + 1 or 2 + 1, therefore 

d T s * (p m dt) n = df * Jlt~ 2 = d s p n e = d t Je + ^ A dt = 0, 
and the continuity equations can be reduced to 

df * (p m dt) n - d T s {*J^) = d, JT^ A dt - d tP n e = 0. 
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So we have 
d t d s D n - d tP n e 

dU T s B n - d T t * (p m dty 



d tP n e - d s {d s H n+l Adt - J e 2 Adt) 



= 



-dj * ( Pm dt) n + d T s {d T s E n+l Adt + *J™ +5 ) 



= 0. 



□ 

Now we show the scheme (10) on the 2D discrete space manifold and the 
time. The second and fourth equations in Eqs.(10) based on Fig.l are 



B n+i _ B n 
At 

->n+l 



+ / 



E n+l _ E n+l 

I * ei | 

rn+l|„ I I rj-re+li 



m+li 



+ J el — 



At 



IP: 



With relations (8), Eqs.(ll) can be rewritten into an implicit scheme (12) 
for TM wave. 



E 1 — E 1 n+± n l \e\\+H 2 |e 2 | + H 3 1 63 1 

£ — + J el = idI 



MO 



At 

trn+l TTn 

-r -J, 



At 



ml 



I * ei I 



► TM (12) 



General schemes 

For real world materials, the constitutive relations are not simple propor- 
tionalities, except approximately. The relations can usually still be written: 

D = eE B = pH, 



but e and p are not, in general, simple constants, but rather functions. With 
Ohm's law , ^ 

E = — J J m = — H, 

a a m 

where a is the electrical conductivity and a m is magnetic conductivity. The 
IDEC schemes can be written as 



- e? e™ +1 + e? m +l - Hl} +1 

e— — — L + a— — L = — — ; — 



At 



* ei I 



H n+1 _ R n H n+1 + R n ^ gn+1 ^ | + ^ | + gn+1 ^ | 



At 



> TE 
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4 Stability, convergence and accuracy 

Now, we analyze the stability for scheme (9). The analysis for scheme (12) 
can be done in the same way. Suppose the fields to be: 

E? +1 = Eft 

H^ 1 = i^cos(fc|* ei |)£, 

where £ is the growth factor of time, and k is the spatial frequency spectrum. 
Substituting (13) into scheme (9), we obtain 



(14) 



= E? + - i ^- ; (l-coa(k\*e 1 \))H?Z 
e\ * e±\ 

H?S = H?--^(E?\ ei \ + E2\e 2 \ + E2\e 3 \){i 

MKll 

Rewrite the first equation of Eqs.(14) as 

e\ * ei|(£ - 1) 

and substitute it into the second equation of Eqs.(14) to obtain 

Hft = H?-^^({l-cos(k\*e 1 \))H?-^ + {l-cos{k\*e 2 \))H{ 
|n | V | * ei| 

+(l-cos(fc|*e 3 |))i?r N ^ ^ 



* e 3 | y £ - 1 

Therefore, we obtain a quadratic equation for £ as follows: 

(l + M )£ 2 -2£ + l = 0, (15) 

where 

M = f(l-cos(fc| *ei|))^i 1 -!^- r + (l-cos(A;| *e 2 \))fa-. 



Pi | V I * e l| I * e 2\ 

+(1 - cos(fe| * e 3 \))(j) 1 -^- J > 

l* e 3|y 



The discriminant of Eq.(15) is 



4-4(1 + Af) = -AM < 0. 

So 

That is to say scheme (9) is unconditional stability. 

By the definition of truncation error, the exact solution of Maxwell's 
equations satisfy the same relation as IDEC scheme except for an additional 
term 0((At) 2 + At\ * e|). This expresses the consistency, and so conver- 
gence for IDEC scheme by Lax equivalence theorem (consistency + stability 
= convergence) . The derivative of IDEC scheme are approximated by first 
order difference. Equivalently, H and E are approximated by linear interpo- 
lation functions. Consulting the definition about accuracy of finite volume 
method, we can also say that IDEC scheme has first order temporal and 
spacial accuracy. 



5 Implementation 

The IDEC scheme of Maxwell's equations was implemented on C++ plat- 
form, consisting of the following steps: 

1. Set the simulation parameters. These are the dimensions of the com- 
putational grid and the size of the time step, etc. 

2. Set the propagating media parameters. 

3. Initialize the mesh indexes. 

4. Assign current transmitted signal. 

5. Compute the value of all spatial nodes and temporarily store the result 
in the circular buffer for further computation. 

6. Visualize the currently computed grid of spatial nodes. 

7. Repeat the whole process from the step 4, until reach the desired time. 
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Fig.3 and Fig.4 exhibit the propagation of Gaussian pulse on rabbit and 
sphere simulated by IDEC. 




Figure 5 
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